Desire lines within each NYC Borough between their NTAs

library(sf)
library(geojsonsf)
library(tidyverse)
library(tmap)

Data pre-processing

Import data from files - NYC NTA Polygons - NYC NTA to tract equivalence - NYS OD

options(scipen = 999)
nyc_nta_geo <- geojson_sf('./data/nyc_nta_2010.geojson')
nyc_tract_nta_equiv <- readxl::read_xlsx('./data/nyc_2010_census_tract_nta_equiv.xlsx')
ny_origin_dest <- read_csv("./data/ny_od_main_JT00_2019.csv")

Define places of interest

boros <- c("Bronx", "Brooklyn", "Manhattan", "Queens")
county_codes <- c("005", "047", "061", "081")
parks <- c("BX10", "BX99", "BK99", "MN99", "QN99")

Create equivalency for counties and tracts with NTAs

county_tract_nta_equiv <- readxl::read_xlsx('./data/nyc_2010_census_tract_nta_equiv.xlsx') %>%
  filter(!(nta_code %in% parks)) %>%
  mutate(county_tract = str_c(`county_code`, `census_tract`)) %>%
  select("county_tract", "nta_code")

Reduce NYS origin destination data to only ntas of interest

nta_ods <- ny_origin_dest %>%
  filter(
    str_sub(as.character(w_geocode), 3, 5) %in% county_codes &
      str_sub(as.character(h_geocode), 3, 5) %in% county_codes
    ) %>%
  mutate(w_county_tract = str_sub(as.character(w_geocode), 3, 11)) %>%
  mutate(h_county_tract = str_sub(as.character(h_geocode), 3, 11)) %>%
  select(h_county_tract, w_county_tract, S000) %>%
  left_join(county_tract_nta_equiv, c("h_county_tract" = "county_tract")) %>%
  rename(h_nta_code = nta_code) %>%
  left_join(county_tract_nta_equiv, c("w_county_tract" = "county_tract")) %>%
  rename(w_nta_code = nta_code) %>%
  mutate(od = str_c(h_nta_code, w_nta_code)) %>%
  group_by(od) %>%
  summarise(
    h_nta_code,
    w_nta_code,
    S000 = sum(S000),
    ) %>%
  unique()

Look exclusively at trips that start and end in the same Borough

intra_borough_nta_ods <- nta_ods %>%
  filter(
    str_sub(h_nta_code, 1, 2) == str_sub(w_nta_code, 1, 2)
    ) %>% 
  mutate(boro = str_sub(h_nta_code, 1, 2))

Reduce NYC geography to NTAs of interest

nta_polys <- nyc_nta_geo %>%
  filter(CountyFIPS %in% county_codes) %>%
  filter(!(NTACode %in% parks)) %>%
  select(BoroName, NTACode, geometry)

Find point on surface

nta_points <- nta_polys %>%
  mutate(geometry = st_point_on_surface(geometry))

## Warning in st_point_on_surface.sfc(geometry): st_point_on_surface may not give
## correct results for longitude/latitude data

Define most popular destinations

intra_borough_dests <- intra_borough_nta_ods %>%
  group_by(w_nta_code) %>%
  summarise(
    w_nta_code,
    S000 = sum(S000)
  ) %>%
  unique() %>%
  left_join(nta_polys, c("w_nta_code"= "NTACode")) %>%
  st_as_sf()

Define most popular origins

intra_borough_origins <- intra_borough_nta_ods %>%
  group_by(h_nta_code) %>%
  summarise(
    h_nta_code,
    S000 = sum(S000)
  ) %>%
  unique() %>%
  left_join(nta_polys, c("h_nta_code" = "NTACode")) %>%
  st_as_sf()

Define desire lines for origins and destinations - Origin and Destination NTA must be different

intra_boro_od_lines <- intra_borough_nta_ods %>%
  filter(h_nta_code != w_nta_code) %>%
  left_join(nta_points, c("h_nta_code" = "NTACode")) %>%
  rename(h_geometry = geometry) %>%
  left_join(nta_points, c("w_nta_code" = "NTACode")) %>%
  rename(w_geometry = geometry) %>%
  mutate(geometry = st_union(h_geometry, w_geometry)) %>%
  mutate(geometry = st_cast(geometry, "LINESTRING")) %>%
  select("od","boro", "S000", "geometry") %>%
  st_as_sf()

Manhattan

Origins

mn_origins <- intra_borough_origins %>%
  filter(BoroName == "Manhattan")

tm_shape(mn_origins) +
  tm_polygons(
    col = "S000",
    style = "jenks",
    title = "Residents"
  ) +
  tm_layout(
    title = "Most popular origins Manhattan",
    legend.outside = TRUE
  )

Desinations

mn_dests <- intra_borough_dests %>%
  filter(BoroName == "Manhattan")

tm_shape(mn_dests) +
  tm_polygons(
    col = "S000",
    style = "jenks",
    title = "Employers"
  ) +
  tm_layout(
    title = "Most popular destinations in Manhattan",
    legend.outside = TRUE
  )

Desire lines

mn_outlines <- nta_polys %>%
  filter(BoroName == "Manhattan")

mn_od_lines <- intra_boro_od_lines %>%
  filter(boro == "MN")

tm_shape(mn_outlines) +
  tm_polygons(
    col = "BoroName",
    title = "Borough",
  ) +
  tm_shape(mn_od_lines) +
    tm_lines(
      col = "#212121",
      lwd = "S000",
      ) +
  tm_layout(
    legend.outside = TRUE,
    legend.text.size = 0.55,
  )

Bronx

Origins

bx_origins <- intra_borough_origins %>%
  filter(BoroName == "Bronx")

tm_shape(bx_origins) +
  tm_polygons(
    col = "S000",
    style = "jenks",
    title = "Residents"
  ) +
  tm_layout(
    title = "Most popular origins in the Bronx",
    legend.outside = TRUE
  )

Desinations

bx_dests <- intra_borough_dests %>%
  filter(BoroName == "Bronx")

tm_shape(bx_dests) +
  tm_polygons(
    col = "S000",
    style = "jenks",
    title = "Employers"
  ) +
  tm_layout(
    title = "Most popular destinations in the Bronx",
    legend.outside = TRUE
  )

Desire lines

bx_outlines <- nta_polys %>%
  filter(BoroName == "Bronx")

bx_od_lines <- intra_boro_od_lines %>%
  filter(boro == "BX")

tm_shape(bx_outlines) +
  tm_polygons(
    col = "BoroName",
    title = "Borough",
  ) +
  tm_shape(bx_od_lines) +
    tm_lines(
      col = "#212121",
      lwd = "S000",
      ) +
  tm_layout(
    legend.outside = TRUE,
    legend.text.size = 0.55,
  )

Brooklyn

Origins

bn_origins <- intra_borough_origins %>%
  filter(BoroName == "Brooklyn")

tm_shape(bn_origins) +
  tm_polygons(
    col = "S000",
    style = "jenks",
    title = "Residents"
  ) +
  tm_layout(
    title = "Most popular origins in Brooklyn",
    legend.outside = TRUE
  )

Desinations

bn_dests <- intra_borough_dests %>%
  filter(BoroName == "Brooklyn")

tm_shape(bn_dests) +
  tm_polygons(
    col = "S000",
    style = "jenks",
    title = "Employers"
  ) +
  tm_layout(
    title = "Most popular destinations in Brooklyn",
    legend.outside = TRUE
  )

Desire lines

bn_outlines <- nta_polys %>%
  filter(BoroName == "Brooklyn")

bn_od_lines <- intra_boro_od_lines %>%
  filter(boro == "BK")

tm_shape(bn_outlines) +
  tm_polygons(
    col = "BoroName",
    title = "Borough",
  ) +
  tm_shape(bn_od_lines) +
    tm_lines(
      col = "#212121",
      lwd = "S000",
      ) +
  tm_layout(
    legend.outside = TRUE,
    legend.text.size = 0.55,
  )

Queens

Origins

qn_origins <- intra_borough_origins %>%
  filter(BoroName == "Queens")

tm_shape(qn_origins) +
  tm_polygons(
    col = "S000",
    style = "jenks",
    title = "Residents"
  ) +
  tm_layout(
    title = "Most popular origins in Queens",
    legend.outside = TRUE
  )

Desinations

qn_dests <- intra_borough_dests %>%
  filter(BoroName == "Queens")

tm_shape(qn_dests) +
  tm_polygons(
    col = "S000",
    style = "jenks",
    title = "Employers"
  ) +
  tm_layout(
    title = "Most popular destinations in Queens",
    legend.outside = TRUE
  )

Desire lines

qn_outlines <- nta_polys %>%
  filter(BoroName == "Queens")

qn_od_lines <- intra_boro_od_lines %>%
  filter(boro == "QN")

tm_shape(qn_outlines) +
  tm_polygons(
    col = "BoroName",
    title = "Borough",
  ) +
  tm_shape(qn_od_lines) +
    tm_lines(
      col = "#212121",
      lwd = "S000",
      ) +
  tm_layout(
    legend.outside = TRUE,
    legend.text.size = 0.55,
  )

## Legend labels were too wide. Therefore, legend.text.size has been set to 0.49. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.